/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/operation.hxx"

#include "linal/device/host/operation.hxx"

namespace linal {

  // Size of T drive the QR_inc
  // The dimension of A and T should be same
  int qr_inc(Hier_ A, Hier_ T) {
    if (!A.get_m() || !A.get_n()) return true;

    /*
    LINAL_ERROR_DETAIL( (A.get_m() == T.get_m()) &&
                        (A.get_n() == T.get_n()) ,
                        "%s\n dim A = ( %4d, %4d )\n dim T = ( %4d, %4d )\n",
                        LINAL_ERROR_MATRIX_DIMENSION_MISMATCH,
                        A.get_m(), A.get_n(), 
                        T.get_m(), T.get_n() );
    */

    FLA_Obj ATL,   ATR,      A00, A01, A02, 
            ABL,   ABR,      A10, A11, A12,
                             A20, A21, A22;

    FLA_Obj TTL,   TTR,      T00, T01, T02, 
            TBL,   TBR,      T10, T11, W12,
                             T20, T21, T22;

    FLA_Part_2x2( ~A, &ATL, &ATR,
                      &ABL, &ABR, 0, 0, FLA_TL );
    FLA_Part_2x2( ~T, &TTL, &TTR,
                      &TBL, &TBR, 0, 0, FLA_TL );

    // factorize based on the dimension of TBR 
    while ( FLA_Obj_min_dim( TBR ) > 0 ){
      FLA_Repart_2x2_to_3x3( ATL,  ATR,       &A00,  &A01, &A02,
                                              &A10,  &A11, &A12,
                             ABL,  ABR,       &A20,  &A21, &A22,
                             1, 1, FLA_BR );
      FLA_Repart_2x2_to_3x3( TTL,  TTR,       &T00,  &T01, &T02,
                                              &T10,  &T11, &W12,
                             TBL,  TBR,       &T20,  &T21, &T22,
                             1, 1, FLA_BR );

      //------------------------------------------------------------
      Hier_ 
        A_11(A11), A_21(A21), A_12(A12), A_22(A22), 
        T_11(T11), T_21(T21), W_12(W12); 

      FLA_Obj 
        &WA_11 = ((FLA_Obj*)A_11.get_buffer())[0],
        &WT_11 = ((FLA_Obj*)T_11.get_buffer())[0];
        
      int b = FLA_Obj_min_dim( WA_11 );

      // dimension of T is modified at this time
      WT_11.m = b;
      WT_11.n = b;
        
      FLA_QR_UT( WA_11, WT_11 );

      apply_q( FLA_LEFT, FLA_CONJ_TRANSPOSE, 
               FLA_FORWARD, FLA_COLUMNWISE,
               A_11, T_11, W_12, A_12);

      if (T_21.get_m() && T_21.get_n()) {
	
        void
          *buf_A_21 = A_21.get_buffer(), *buf_A_12 = A_12.get_buffer(),
          *buf_A_22 = A_22.get_buffer(), *buf_W_12 = W_12.get_buffer(),
          *buf_T_21 = T_21.get_buffer();
	
        int
          cs_A_21 = A_21.get_cs(), cs_A_12 = A_12.get_cs(),
          cs_A_22 = A_22.get_cs(), cs_W_12 = W_12.get_cs(),
          cs_T_21 = T_21.get_cs();
        
        for (int i=0;i<T_21.get_m();++i) {

          // make T_21 square
          LINAL_ELEM_FLA(_T_21, i,0).m = b;
          LINAL_ELEM_FLA(_T_21, i,0).n = b;

          FLA_QR2_UT_internal( WA_11,
                               LINAL_ELEM_FLA(_A_21, i,0), 
                               LINAL_ELEM_FLA(_T_21, i,0),
                               fla_qr2ut_cntl_leaf );

          for (int j=0;j<W_12.get_n();++j) {

#pragma omp task firstprivate(A_21, T_21, W_12, A_12, A_22, i,j)
            {
              FLA_Apply_Q2_UT_internal( FLA_LEFT, FLA_CONJ_TRANSPOSE,
                                        FLA_FORWARD, FLA_COLUMNWISE,
                                        LINAL_ELEM_FLA(_A_21, i,0), 
                                        LINAL_ELEM_FLA(_T_21, i,0),
                                        LINAL_ELEM_FLA(_W_12, 0,j), 
                                        LINAL_ELEM_FLA(_A_12, 0,j),
                                        LINAL_ELEM_FLA(_A_22, i,j),
                                        fla_apq2ut_cntl_leaf );
            }
          }

#pragma omp taskwait

        }
      }
	
      //------------------------------------------------------------
      FLA_Cont_with_3x3_to_2x2( &ATL,  &ATR,       A00, A01,  A02,
                                                   A10, A11,  A12,
                                &ABL,  &ABR,       A20, A21,  A22,
                                FLA_TL );
      FLA_Cont_with_3x3_to_2x2( &TTL,  &TTR,       T00, T01,  T02,
                                                   T10, T11,  W12,
                                &TBL,  &TBR,       T20, T21,  T22,
                                FLA_TL );
    }

    return true;
  }
}
